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1997|) . This layer, called tachocline, is supposed to play 



Inferring the equatorial solar tachocline from frequency 
splittings 

T. Corbard, G. Berthomieu, J. Provost, and P. Morel 

Laboratoire G.-D. Cassini, CNRS UMR 6529, Observatoire de la Cote d'Azur, BP 4229, 06304 Nice Cedex 4, FRANCE 
Received 18 July 1997; accepted 24 October 1997 

Abstract. Helioseismic inversions, carried out for several 
years on various ground-based and spatial observations, 
have shown that the solar rotation rate presents two prin- 
cipal regimes: a quasi-rigid rotation in the radiative inte- 
rior and a latitude-dependent rotation in the whole con- 
vection zone. The thin layer, named solar tachocline, be- 
tween these two regimes is difficult to infer through in- 
verse techniques because of the ill-posed nature of the 
problem that requires regularization techniques which, in 
their global form, tend to smooth out any high gradient 
in the solution. Thus, most of the previous attempts to 
study the rotation profile of the solar tachocline have been 
carried out through forward modeling. In this work we 
show that some appropriate inverse techniques can also 
be used and we compare the ability of three ID inverse 
techniques combined with two automatic strategies for 
the choice of the regularization parameter, to infer the 
solar tachocline profile in the equatorial plane. Our work, 
applied on LOWL (LOWL is an abbreviation for low de- 
gree denoted by L) two years dataset, argue in favor of 
a very sharp (0.05 ± 0.03i?©) transition zone located at 
0.695 ± 0.005-R© which is in good agreement with the 
previous forward analysis carried out on Global Oscilla- 
tions Network Group (GONG), Big Bear Solar Observa- 
tory (BBSO) and LOWL datasets. 

Key words: Sun: interior - Sun: oscillations - Sun: rota- 
tion methods: numerical 



1. Introduction 

Helioseismic inversions of the solar p-modes frequencies 
splitted by rotation have shown that there is, at the base 
of the convection zone, a thin transition layer separating 
two regimes of rotation, a strong differential rotation in 
the convection zone and a quasi rigid rotation in the ra- 
diative interior (e.g. Thompson et al. 1996; Corbard et al. 



an important role in the solar dynamo, in the transport 
of angular momentum and in the mixing of chemical ele- 
ments. Its position r c and thickness w give constraints to 
the theo ries d escribing its str ucture and evolution (Spiegel 
& Zahn |1992| ; Gough & Sekii |1997| ). Different estimations 
of these parameters have been obtained so far mostly by 
using forw ard m ethod s (Kosovichev 1996 ; Charbonneau 
et al. |!997t Basu |l997| ). 

The aim of this work is to test and compare the abil- 
ity of some inversion methods to infer the location and 
the width of the solar tachocline, and then to apply these 
methods to helioseismic data. We compare three ID least- 
squares methods. They differ essentially by the mean used 
to regularize the ill-posed inverse problem of inferring the 
equatorial solar rotation rate from the observed frequency 
splittings. The first method is the most commonly used 
Regularized Least-Squares (RLS) method with Tikhonov 
regularization ( Tikhonov fc Arscnin 197*^ ), the second one 
is the Modified Truncated Singular Value Decomposi- 
tion (MTSVD) introduced by Sekii and Shibahashi fll988| ) 
which uses a regularization term of the same form but with 
a discrete truncation parameter instead of the continuous 
Tikhonov regularization parameter. The third method, in- 
troduced by Hansen & Mosegaard (1996), is called Piece- 
wise Polynomials TSVD (PP-TSVD) and is a modification 
of the MTSVD method that can preserve discontinuities 
of the solution. 

In Sect. [2J we briefly recall the inverse problem and 
define our parameterization of the tachocline. Section]^ 
gives the two strategies studied in this work for inferring 
the rapid variation of the rotation. We test these methods 
by inverting artificial data in Sect. ^ and then, in Sect. ^, 
we use this study in order to infer the location and thick- 
ness of the solar tachocline in the equatorial plane from 
data observed by the LOWL instrument (Tomczyk et al. 



1995) 
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2. Direct analysis and parameterization of the 
tachocline 

Frequency splittings Av nim = v nim - i/ n ;-m between 
modes with the same radial order n and degree I but differ- 
ent azimuthal orders m are induced by the solar rotation 
f2(r, 9) expressed as a function of the radius r and colat- 
itude 9. For a slow rotation, assumed to be symmetric 
about the equator, and moderate or high degree modes, 
these splittings are given by: 

Av nlm =m [ 2 [ °K nl {r)Pl n {cos9) 2 fl(r,6) sin dr <ffl,(l) 
Jo Jo 

where K n i (r) are the so-called rotational kernels that can 
be calcu lated for each mode from a solar model (Morel et 
al. 1997). In the following, they are assumed to be known 
exactly. There exists additional terms that are not taken 
into account in Eq. ([!]) but, as discussed in Corbard et al. 
( |1997| ), they do not influence inversion above 0.4R Q . As 
the aim of this work is not to sound the rotation of the 
core, Eq. ([!]) is a good approximation. P[ n (cos9) are nor- 
malized Legendre functions. Their asymptotic property 
leads, as discussed by Antia et al. ( 1996| ), to the follow- 
ing expression that shows the sectoral (i.e. I — m) modes 
splittings as weighted averages of the equatorial rotation 
rate Q eq {r) = n(r,90°): 



I / K nl (r) n eq (r) dr. 



(2) 



We note that the validity of this ID approximation is l- 
dependent. Indeed, the higher the degree, the more the 
latitudinal kernel P/(cos 9) 2 sin is peaked at the equator. 

Following Charbonneau et al. ( 1997 ), we define the 
location and the width of the transition zone in the equa- 
torial plane as the parameters r c and w respectively of the 
following erf function which fits the rotation law in this 
plane: 



n eq (r) =Cl + - fio) ( 1 + erf (^—^ 
2 V V u.bw 



(3) 



Here Clo and Cti represent the mean values of the rota- 
tion in the radiative interior and in the convection zone 
respectively. 

In order to compare different ID inverse methods, we 
have built several sets of theoretical sectoral frequency 
splittings that correspond to different given rotation laws 
with fixed parameters r c , w, £Iq, £li but with a function of 
the colatitude in order to mimic the latitudinal differential 
rotation of the convection zone: 



n(r, 9)=no+-(n 1 ~Acos 2 9-Bcos 4 9-n )(l+ erf(- — - 
2 V \0.5ui 



(4) 



w = w, CIq = Qq and Cl\ = Q-i . We compute the splittings 
Ap n u from Eq. (Q) for a set of modes corresponding to the 
set of LOWL data used in Corbard et al. (1997) and we 
add a normally distributed noise Sv n u G Af(0,a„i). For 
each mode (n, I) the standard deviation of the noise a n i 
has been taken equal to: 



(Tnl 



(5) 



Evidently, for any choice of constants A and B, the 
searched parameters for these rotation laws are f c = r c , 



where a n i is the error derived from the observers' uncer- 
tainties for a splitting Av> n u, and k a is an integer used to 
vary the level of the noise that we introduce in the data. 
Doing this, we take into account the fact that the error ob- 
tained on the observed splitting varies with the frequency 
and the degree of the mode which is certainly more real- 
istic than taking the same average standard deviation for 
all the modes. From those noisy splittings, the equatorial 
rotation profile is obtained by inverting Eq. (E) and this 
profile is then fitted by the erf function Eq. (3) leading 
to the parameters f c , w, n , £l\ which will be compared 
to the initial parameters. 



3. Strategies for inferring rapid variations of the 
rotation 

The three inverse methods used in this work are detailed 
in Appendix |A| They all use a grid of 50 points in radius 
distributed according to the density of turning points of 
observed modes. The most important difficulty in infer- 
ring the thickness of the tachocline from inverse methods 
results from the fact that the problem of solving Eq. (^) 
is an ill-posed problem and this is strengthened by the 
fact that rotational kernels give redundant information 
about the outer part of the sun whereas they have only 
low amplitude in the solar core for the observed mode set. 
Numerically, this produces a high value for the condition 
number (defined as the maximum singular value divided 
by the smallest singular value) of the discretized problem 
Eq. ( |A5| ) (typically A max / A min ~ 2xl0 8 in our implemen- 
tation) and the singular values decay rapidly. This high 
value of the condition number means that the solution 
of the initial problem is highly sensitive to the numerical 
errors and the noise contained in the data. Therefore we 
have to introduce some a-priori knowledge on the rotation 
profile. Unfortunately this regularization tends to smooth 
out every rapid variation in the solution. By using global 
regularization, we make the implicit assumption that the 
real rotation is smooth everywhere and therefore the in- 
formation about the thickness of a rapid variation of the 
rotation profile is not directly readable from the solutions 
obtained by classic inversions. There are however several 
ways for overcoming these difficulties. 
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3.1. Local deconvolution of the result obtained from linear 
inversions: the use of averaging kernels 

The first way is to have a good understanding of the pro- 
cess by which the inversion smoothes the solution: using 
this information, we may be able to inverse this process 
and to acquire a more realistic view of the rotation. This is 
what Charbonneau et al. ( 1997] ) have done in combination 
with the so-called Subtractive Optimal Localized Average 
(SOLA) (Pijpers & Thompson |1992] , |1994[ ) method. This 
can be generalized for any linear inversion as RLS method 
used in this work. The solution fi(ro) obtained at a tar- 
get location ro can be viewed as a weighted average of 
the 'true rotation' f2(r), the weighting function being the 
averaging kernel fc(r, ro) that can always be estimated at 
any r : 



Appendix Al. We estimate that this cannot be made for 



fi(r ) = / «(r, r )fi(r) dr. 
Jo 



(6) 



If we suppose that the averaging kernels obtained at any 
depth can be approximated by a translation of the aver- 
aging kernel obtained at the middle of the transition i.e. 
n(r, f c ), then we can define k c by n c (r — r c ) = n(r, f c ) and 
Eq. (||) reduces to a convolution equation: 



O(ro) = / K c(r - ro)O(r) dr ^> Q(r) = n c {r) * tt(r) (7) 
Jo 

Finally, if the 'true rotation' can be well approximated by 
an erf function of the form given by Eq. (^), and if we 
approximate the kernel n c (r — ro) by a Gaussian function 
of the form: 



K c (r - r ) ^ exp \-(r - r ) 2 /A 2 r 



(8) 



then the inferred solution is also an erf function of the 
form Eq. (§ but with a larger width w. A simple decon- 
volution gives the following relation between the searched 
width w and the inferred width w: 



= w c = y/w 2 - 4A2, 



(9) 



which defined the corrected inferred width w c . 

This result is valid only under a large number of as- 
sumptions that may be quite distant from the reality. Es- 
pecially the reduction to a convolution form is certainly 
not valid because of the extent of averaging kernels that 
tend to increase rapidly toward the solar core. Moreover 
the profile of the rotation rate may be much more compli- 
cated than a simple erf function. However, the tachocline 
is thin and the averaging kernels have nearly the same 
profile in its whole extent. Thus this is certainly a good 
approach to get a quantitative idea of how the inversion 
enlarges the 'true rotation' transition. We note that if we 
obtain A r > w/2 this certainly means that some of the 
previous assumptions are not valid. In this work, we have 
applied this 'deconvolution method' on the solutions ob- 
tained by Tikhonov inversions computed as explained in 



MTSVD method because the corresponding averaging ker- 
nels are less well peaked and exhibit a more oscillatory 
behavior (see Fig. ^ hereafter). 

3.2. Non linear regularization 

The second way to estimate the location and thickness 
of the tachocline, is to build inverse methods that are 
capable of producing solutions with steep gradients. The 
idea is to apply a local regularization instead of the global 
Tikhonov regularization term. This leads to a non linear 
problem and piecewise smooth solutions. This approach 
has recently found useful applications in image process- 
ing for edge-preserving regularization (Aubert et al. 
and total variation (TV) denoising (Vogel & Oman 



1994) 



1996 



1997 ). In particular, the TV of / is defined as the 1-norm 



of the first derivative of / and this is the definition of 
smoothness that we use in the PP-TSVD inverse method. 
Therefore, the results obtained by this method, detailed in 
Appendix [A.2| , represent a first attempt to use this class 
of inversion with non linear regularization on helioseismic 
data. 

4. Tests with artificial data: results and discussion 

4-.1. The key: how to choose regularization parameters 

Whichever regularized inverse method we use, a very 
important point is the choice of the regularization pa- 
rameter which can be a discrete truncation parameter k 
(MTSVD, PP-TSVD, Eq. flAl^) ) or a continuous parame- 
ter A (Tikhonov, Eq. (|A9|)). This choice is specially impor- 
tant if we want to infer a quantity like the width of a zone 
with high gradients which is directly affected by the reg- 
ularization. Several methods for choosing the regulariza- 
tion parameter have been proposed that tend to establish 
a balance between the propagation of input errors and the 
regularization (see e.g. Badeva & Morozov (1991), Thomp- 
son & Craig ( [1992D and Han sen ( |1992| |1994[ ) for a general 
review and Thompson ( 1992|), B arett (1993) and Stepanov 
& Christensen-Dalsgaard ( 1996| ) for applications in helio- 
seismic inversions). In this work we test and compare the 
ability of two of these automatic strategies, namely the L- 
curve criterion (Hansen |1992| ) and the Generalized Cross 
Validation (GCV) criterion (Wahba |l977|; Golub et al. 



1979), to reproduce a good estimation of the tachocline 



profile from noisy data. 

The importance of the choice of the regularization pa- 
rameter can be illustrated by the following figures (Figs. 
U, ||, ||, |J ) where the results of the fit of the solution by an 
erf function are plotted as a function of the regularization 
parameter. 

Figure |l| represents the variation of the four erf- 
parameters f2o, fii, f c and Q deduced from a Tikhonov 
inversion as a function of the logarithm of the regulariza- 
tion parameter. The four initial parameters were Qo = 425 
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Log 10 (\) Log, (A) 

Fig. la-d. Inferred parameters Clo, Cli, f c , w and corrected 
inferred parameter w c against the logarithm of the Tikhonov 
regularization parameter A. Error bars result from the fit of 
the solution by an erf function taking into account the propa- 
gation of noise through the inverse process but not the existing 
correlations between the results obtained at two different ra- 
dius. The initial parameters are indicated by dashed lines. The 
GCV and L-curve choices are shown by the full star and the 
circle respectively. The input rotation law was not dependent 
on the latitude (A = B — 0) and the level of noise was small 
(ka = 10). 



convection zone close to that derived by 2D inversions. We 
have set A = 55 nHz and B = 75 nHz which are mean 
values derived from observations of the plasma motion at 
the solar surface (Snodgrass & Ulrich |l990| ) . This choice 
for the input rotation law and errors is referred as the 
'realistic case' in the following. The Eq. (|l|) with m = I has 
been used to compute the frequency splittings of sectoral 
modes and ID Tikhonov inversions have been performed 
again in order to infer the equatorial rotation rate from 
Eq. (§. 




-6 -4 -2 -6 -4 -2 

Log 10 (\) Log, (A) 



nHz, fli = 460 nHz, r c = O.69i? and w = 0.05i? Q . 
In this case, called the 'ideal case' in the following, the 
added errors were small (£^=10) and the initial rotation 
law was not dependent on the latitude {A = B = 0). 
The choices designated by L-curve and GCV strategies 
are shown by the full star and the circle respectively. In 
addition we have plotted the corrected inferred width w c 
given by Eq. @ and computed by calculating systemati- 
cally the averaging kernel at ro = f c (as shown on Fig. ||a 
for the GCV choice). The GCV criterion leads always to 
a lower regularization than the L-curve choice and then 
tends to reduce the smoothing of the solution. In most of 
our tests, as in Figs. |l|a, c, d, the GCV choice corresponds 
to a point where the errors deduced from the fit become 
small whereas the L-curve criterion gives a point beyond 
which a rapid variation of the fitted parameters with in- 
creasing regularization occurs. The fact that the values of 
the fitted parameters are nearly constant between these 
two points shows that, for this level of noise, the method 
is robust in that sense that the choice of the precise value 
of the regularization parameter is not a crucial point: any 
choice that tends to establish a balance between the prop- 
agation of input errors and the regularization is able to 
produce good results. 

Let us now look at the behavior of this method for a 
more realistic example. For this we take a level of noise 
similar to the one given by observers (k a = 1) and we 
build frequency splittings of sectoral modes by taking into 
account a latitudinal variation of the rotation rate in the 



Fig. 2a-d. The same as in Fig |]j but with more realistic input 
errors (fc CT = 1) and an input rotation profile with latitudinal 
variation in the convection zone (A = 55 nHz, B — 75 nHz). 



Figure y represents the results of these inversions in 
the same form as Fig. Q and for the same initial erf- 
parameters. There are two essential points to be seen on 
this figure. The parameter flo in Fig. ||a is systematically 
under-estimated of about 4 nHz. A detailed analysis shows 
that this effect is strongly related to the introduction of 
a latitudinal variation of the rotation rate in the convec- 
tion zone. The assumption, used in the ID inversions, that 
sectoral modes are sensitive only to the equatorial com- 
ponent of the rotation rate is not valid for low degree I 
modes (e.g. Antia et al. 1996) , and these modes sound 



the deep interior. This may explain some perturbation for 
the determination of the parameter Qq that represents the 
mean value of the rotation rate in the radiative interior. 
The difference between splittings of sectoral modes com- 
puted from Eq. (^) and Eq. (|l|) is below 1 nHz for the 
observed modes having their turning points above 0.4_Rq. 
The large resulting difference in £Iq is due to the fact that 
high I sectoral modes see only the equatorial rotation rate 
and then fix the inferred value fir equal (or nearly equal 
as in Fig. ||b) to the initial value Oi while lower degrees 
sectoral modes are sensitive to the differential rotation of 
the convection zone and this effect can only be accounted 
for in the inverse rotation law by a substantial lowering in 
f2o- Furthermore wc have checked that two rotation laws 



T. Corbard et al.: Inferring the equatorial solar tachocline from frequency splittings 



5 



with the same fii but with a difference of 4 riHz in f2o and 
two rotation laws with the same f2o but with or without 
latitudinal variation in the convection zone, induce a dif- 
ference of the same order in the sectoral modes frequency 
splittings. 

The second important point is that, in Figs, ^c, d, 
the estimation w of the width of the tachocline increases 
rapidly between the GCV and the L-curve points whereas 
its location f c decreases rapidly from 0.688J?© down to 
0.674i? Q As in Fig. [l]d, the deconvolution made by us- 
ing averaging kernels tends to correct this behavior for 
the estimation of the width but, in this case, the GCV 
choice remains over-estimated for about O.O15i?0 and the 
L-curve choice is still very distant from the initial value. 
Tests made with different input parameters show that, as 
in Figs. ||c, d and for that level of noise, the GCV choice is 
always better than the L-curve choice for the estimation 
of the location and the width of the tachocline. This point 
will be illustrated and discussed in the next section for the 
estimation of widths between 0.03 and O.lli?©. 




Fig. 3a-d. The same as in Fig. ^ ('ideal case') but for MTSVD 
(full line) and PP-TSVD (dashed line) methods and against 
the truncation parameter k. The L-curve choice for MTSVD 
method is outside the plot on panel b. 

Similar figures (Figs. |, |) can be plotted for MTSVD 
and PP-TSVD methods where the continuous regulariza- 
tion parameter is replaced by the discrete truncation pa- 
rameter. Results obtained in the 'realistic case' (Fig. [|) 
have again a larger dispersion and exhibit the same sys- 
tematic deviation for the determination of fio- Another 
interesting point is that, as shown on Figs. |3|d, |]d and also 
in the next section, the PP-TSVD method tends to give 
an under-estimation of the width whereas the MTSVD 
method tends to give an over-estimation of this parame- 
ter. This may be very useful in order to give a bounded 
estimation of the true width. For these two methods, the 
choice of the optimal truncation parameter k through the 
L-curve criterion needs the evaluation of the curvature of 




10 20 30 10 20 30 

k k 



Fig. 4a-d. The same as in Fig. g ('realistic case') but for 
MTSVD (full line) and PP-TSVD (dashed line) methods and 
against the truncation parameter k. The L-curve choice for 
MTSVD method is outside the plot on panels b,c and d. 

discrete L-curve. This can be done carefully by an appro- 
priate 2D curve fitting. Nevertheless our experience shows 
that it is difficult to do this systematically with the same 
fit procedure for any level of noise and input rotation law. 
Furthermore, when this is done carefully, this choice leads 
to results for the tachocline profile that are always worse 
than the ones obtained from the GCV choice. Thus, in the 
following, results are shown only with the GCV criterion 
for MTSVD and PP-TSVD methods. 



™=0.05 R , r c =0.69 R G 
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Fig. 5a-c. Solutions obtained between 0.4 and 0.8-Rq from the 
three inverse methods with the GCV choice of regularization 
parameters. The input rotation law was the same as in Figs. ^| 
('realistic case'). The equatorial component of the initial law 
is shown by dashed line whereas the fits of the inverse solutions 
are shown by full lines. 

Figure || shows the solutions obtained from the three 
methods with the GCV choices indicated on Figs. | and 
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m. The error bars on the PP-TSVD method (Fig. |c) were 
obtained by assuming that the method is linear i.e. the 
dependence of H (defined in Eq. ( Al6[ )) relatively to the 
data vector W is neglected. This is indeed not the case 
and a Monte-Carlo approach for estimating errors may be 
more realistic. We note however that the two other meth- 
ods (Tikhonov and MTSVD) are linear only for a given 
regularization parameter. Since this parameter is chosen 
through automatic strategies, it depends also on the data. 
Thus, strictly speaking, these methods are also non-linear 
methods. Nevertheless, the automatic choices are built so 
that they are not too much sensitive to little change in the 
data and that justify the linear approximation. The cor- 



k (r.r ) 




Fig. 6a and b. Averaging kernels computed at ro = f c . For 
Tikhonov method the dashed line represents the Gaussian ap- 
proximation of the kernel used for the local deconvolution of 
the solution shown on Fig. Ha. 
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Fig. 7a-c. Difference between the inferred width and the ini- 
tial width (Sw = w~w) against the initial width for PP-TSVD 
(triangles) and MTSVD (circles) methods, both computed with 
the GCV choice for the truncation parameter. Squares are 
for the Tikhonov method with GCV criterion (full line) and 
L-curve criterion (dashed line). For this latter method we plot 
the difference between the corrected inferred width and the 
initial width (Sw = w c — w). a k a = 10, A = B — as in Figs. 
§, | ('ideal case'); b fc CT = 1, A = B = 0; c fe CT = 1, A = 55, 
B = 75 as in Figs. 0, [| ('realistic case') 



responding averaging kernels computed at r = f c (Fig. 0) 
show that whereas the Gaussian approximation is rather 
good for the Tikhonov method, the large oscillations in the 
convection zone obtained for the MTSVD method make 
difficult the use of a local deconvolution in that case. 

4-2. Tests for width between 0.03 and 0.11 Rq 

An important point is to test the ability of a method to 
give a good estimation of the er /-parameters for a large 
domain of variation of the width of the tachocline. We 
first study in Fig. the behavior of the different methods 
and automatic strategies between the 'ideal case' and the 
'realistic case' for one realization of input errors. Then, in 
Fig. 0, we have carried out a Monte-Carlo approach in or- 
der to have a better estimation of the errors on the widths 
deduced from the fit of the solutions for the 'realistic case'. 

Figure shows the inferred width w (for MTSVD and 
PP-TSVD methods) and the corrected inferred width w c 
(for the Tikhonov method) as functions of the initial width 
w and for one realization of the input errors. Figure 0a 
represents the same example as Figs. 0, || ('ideal case' ), in 
Fig. Qb we increase the level of noise (k a — 1), and finally 



we set an input rotation law with a latitudinal dependence 
in the convection zone so that the Fig. 0c is for the same 
example as Figs. ||, || ('realistic case'). 

In Fig. 0a , the results for w fit the real value within 
0.02^0 except for PP-TSVD and widths above O.9i? , 
and the two regularization procedures (L-curve and GCV) 
give almost the same result. 

The comparison of Figs. 0a and 0b clearly indicates 
that the results obtained for Tikhonov method with the 
L-curve criterion (dashed curves) are very sensitive to the 
level of noise and are not adapted to the actual errors of 
observed data. The deconvolution method using Tikhonov 
inversion with GCV criterion appears to be the less sen- 
sitive to the noise level and the most stable for widths 
between 0.03 and 0.11i? Q . We see again that the results 
obtained from MTSVD and PP-TSVD lead respectively 
to an over-estimation and an under-estimation of the real 
width. Figure 0c illustrates the effect of a latitudinal de- 
pendence of the rotation in the convection zone: an in- 
creasing over-estimation of w from the Tikhonov method 
with GCV criterion and a general larger dispersion of the 
results. 
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existing correlations between the inferred rotation values 
obtained at two different radius are not taken into account 
in the fit of the solution by an er/-function. Secondly, for 
the PP-TSVD method, the non-linearity of the method 
is not taken into account in the estimation of the prop- 
agation of noise through the inverse process. Making the 
fit in the right way, i.e. taking into account correlations, 
may lead to a lower dispersion of the results and then 
our estimation of the error on the inferred widths may be 
over-estimated. Nevertheless, the effects of these two ap- 
proximations are not easy to estimate a priori and need a 
more complete analysis in future work. 

5. Results for LOWL data 



Fig. 8. The same as in Fig. (?]c ('realistic case') but each points 
is the mean value of the results obtained for 500 realizations of 
input errors. Error bars represent a 68.3% confidence interval 
on w. 



In Fig. ||, we have performed 500 realizations of input 
errors for each initial width and each point shown in this 
figure represents the mean value of the 500 inferred or cor- 
rected inferred widths for a given initial width and a given 
method. Error bars represent a 68.3% confidence interval 
which contains the nearest 341 inferred widths from the 
mean value but they are not necessarily symmetric around 
this value. This study shows that the Tikhonov and PP- 
TSVD methods with the GCV criterion are the most re- 
liable for estimating the width in the most realistic case. 
They lead, respectively to an over-estimation and under- 
estimation of the width of about O.Oli?© at the maximum 
for initial widths between 0.03i?© and O.lli?©. In that 
range, the standard deviation obtained for 500 realizations 
of input errors is around 0.02i?© for Tikhonov method 
and much larger (up to 0.05-R© for w = O.lli?©) for PP- 
TSVD method which then appears to be well adapted 
only to infer very sharp transitions. Let u>i represent the 
widths deduced from N r hypothetical (non-observed) re- 
alizations of the unknown true width Co. In the Monte- 
Carlo method we suppose that we can approximate the 
distribution of (Co — u>i, i = 1, ..N r ) by the distribution of 
(uj — uji, i=l, N r ) where oo is the width deduced from 
the observed dataset and u>i are the widths deduced from 
datasets built by setting Co = u in the model. As we can 
not insure that lo is very close to Co, the underlying as- 
sumption is that, in the range of uncertainty concerning 
Co (say 0.03 — O.lli?©), the way in which errors propagate 
through the i nvers e process does not vary rapidly (see e.g. 
Press et al. 1992 ). The fact that, in Fig. |[ error bars 
grow rapidly with the initial width for PP-TSVD method 
makes difficult the use of the Monte-Carlo results for esti- 
mating the statistical behavior of this method. There are 
nevertheless two factors that may introduce bias in these 
estimations of the errors on the inferred widths. First, the 



0.6 



Fig. 9. Equatorial tachocline profiles obtained from LOWL 
data by PP-TSVD (triangles) and Tikhonov (squares) meth- 
ods with GCV criterion. Error bars represent the la errors 
estimated on the solution by assuming the linearity of the in- 
versions. The full and dashed curves represent respectively the 
fit of the PP-TSVD and Tikhonov solutions by an er /-function 
between 0.4 and 0.8-Rq. 



This section gives the results obtained from the two years 
(2/26/94-2/25/96) observati ons by the LOWL instr umen t 
in Hawaii (Tomczyk et al 



1995; Corbard et al. 1997) 



These data contain 1102 modes with degrees up to I — 99 
and frequencies between 1200 and 3500 ^Hz. For each 
mode (n,l), individual splittings are given by, at best, 
five a-cocfficicnts of their expansion on orthogonal poly- 



nomials defined by Schou et al. (1994). For this work, we 
assume that the previous simulations provide an estima- 
tion of the bias introduced by the methods and we use 
these values in order to correct the inferred tachocline pa- 
rameters. This supposes the closeness of the model used 
in the simulation to the reality and a good estimation of 
the errors in the data. Furthermore, we use the sum of 
odd a-coefficients as a first approximation for the sectoral 



splittings i.e. Av n u 



ii i 'ni 
n + a 3 



. This approximation is 



exact for all the rotation laws such that 



l 2j + l 



V j > 2 
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Table 1. Inferred er /-parameters obtained from LOWL data. The L-curve criterion has not been used for methods with discrete 
truncation parameters. 



Methods 


€l (nHz) 




Hi (nHz) 




fc/Re 


W {c) /Rq 




GCV 


L-curve 


GCV 


L-curve 


GCV 


GCV 


Tikhonov 

MTSVD 

PP-TSVD 


429.3 ±0.5 

429.4 ±0.7 
429.6 ± 0.2 


427.9 ±0.3 


457.7 ± 0.3 
457.0 ± 0.5 
456.4 ± 0.3 


460.4 ± 0.4 


0.693 ±0.002 
0.693 ±0.003 
0.693 ±0.009 


0.067 ±0.010 
0.062 ±0.009 
0.031 ±0.017 



(which is the case for the rotation laws Eq. ^ used in our 
model). When this is not the case the latitudinal kernel 
associated to a"' + a% 1 + aV; 1 is less peaked at the equa- 
tor than the one associated to the sectoral splittings (i.e. 
P/(cos#) 2 sin#, see Sect. ||) and thus &i represents a lati- 
tudinal average of the rotation in a larger domain around 
the equator. However the kernel associated to the sum of 
three a-coefficients is less /-dependent. 

Results obtained by the three methods are summarized 
in Table 1. They are in very good agreement for the loca- 
tion of the tachocline and the mean values of the rotation 
rate in the radiative interior and convection zone but more 
dispersive concerning the determination of the width. The 
tests discussed above have shown that this may be related 
to the level of noise contained in the data. The equatorial 
tachocline profiles obtained by Tikhonov and PP-TSVD 
methods with GCV criterion are shown in Fig. [| Accord- 
ing to the previous sections, we use the GCV choice in 
order to infer the location and the width of the equato- 
rial tachocline. Nevertheless, for Slo and fii the L-curve 
choices may be useful in order to see the amplitude of the 
variation of the inferred parameters against the regulariza- 
tion parameter. The errors cited in this Table are just the 
result of the fit of the solution by the er/-function. The 
variation of the inferred erf parameters against the regu- 
larization, as shown by Fig. [lO] for the Tikhonov method, 
and the previous Monte-Carlo simulations can help us to 
estimate error bars that may be more realistic. 

Figure P~0] el shows that the evaluation of the mean value 
of the rotation rate in the radiative interior (^o) is not 
much sensitive to the regularization. Nevertheless, we have 



shown in Sect. 4.1 that this parameter tends to be sys- 
tematically under-estimated of about 4 nHz because of 
the influence of the latitudinal variation of the rotation in 
the convection zone on the low I sectoral splittings. For 



the sum a" 



■ni 



-,nl 



-,nl 



the latitudinal kernel is less l- 



dependent so that this systematic offset may be smaller 
than 4 nHz. We take this effect into account by increas- 
ing the estimation of the error and our final interval for 
this parameter becomes: 427.5 < Clo < 434.5 nHz. The 
mean value of the equatorial rotation rate in the convec- 
tion zone is less subject to systematic errors but may be 
under-estimated by the GCV choice (c.f. Figs. ||b, |Jb, |J). 




Log 1D (A) 



Fig. 10. Variation of the inferred parameters Qo,Qi,f c ,w and 
w c as a function of the logarithm of the regularization parame- 
ter for the Tikhonov inversion of LOWL data. Graph markers 
have the same meaning as in Fig. |l| The L-curve choice of w 
is outside the plot on panel d. 



The difference between the GCV choice and the L-curve 
choice is about 3 nHz on Fig. Thus we estimate that 
f2i = 459.0 ±1.5 nHz. We note that we do not attempt 
to use the points of the solution found under 0.4i?o or 
above O.8-R0 (c.f. Fig. [)]). Therefore &i does not take into 
account the eventual rapid variation of the rotation near 



the surface or at 0.9i?o (Antia et al. 1996) and is not 
sensitive to the core rotation. The ratio q = Q Q /Cli ob- 
tained from helioseismic data is an important test for the 
theories of the tachocline dynamics. Spiegel and Zahn's 
( [1992D theory leads to q = 0.90 whereas Cough's ( [1985D 
one leads to q = 0.96. Our results give 0.93 < q < 0.95 
which is intermediate between the two theoretical esti- 
mates. Similar results have already been pointed out by 



Gough & Sekii ( |1997| ). 

For the estimation of f c , we find in Fig. |l0| that the 
L-curve criterion leads to a lower value than the GCV cri- 
terion as we had found in Fig. ^. As discussed in Sect. LI, 
we think that the GCV choice is more reliable but may 
lead to an under-estimation of about 0.002i? Q . Therefore 
our final estimation for the location of the center of the 
tachocline in the equatorial plane is: r c = O.695±O.OO5J?0. 
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This value, estimated in the equatorial plane, is inter- 
mediate between the two values previously obtained by 
forward methods (c.f. Table 2). We note however that 
whereas our work just look for the equatorial component 
of the tachocline, the previous works assume that the solar 
tachocline presents the same profile at any latitude. This 
may lead to bias if, as suggested by Charbonneau et al. 
( 1997 ) from LOWL data, the tachocline is prolate i.e. is 
located deeper at the equator than at higher latitudes. 

The tests discussed in the previous sections show that 
the L-curve choice is not reliable for the estimation of the 
width and suggest three ways for estimating the width of 
the tachocline from GCV criterion: 

First, the true value is supposed to lie between 
the MTSVD and PP-TSVD estimations. That gives 
0.031i? Q < w < 0.062i? Q . 

Secondly, for the Tikhonov method, since the error 
bars have roughly of the same amplitude in the whole 
range 0.03 — O.lli?© of initial widths (Fig. ||) , we can use 
the Monte-Carlo simulation. Near w = 0.07i?o (the in- 
ferred value reported in Table 1 being w = 0.067), Fig. [s] 
shows that the Tikhonov method leads in mean to a sys- 
tematic over-estimation of about 0.005-R© with a 68.3% 
confidence interval around ±0.02i? Q . Thus we obtain by 
this way w ~ 0.062 ± 0.020i?©. 

Thirdly, the PP-TSVD method is though to pro- 
duce, in mean, an under-estimation of the width of about 
O.Oli?© but with a larger dispersion of the results for the 
large widths so that we are not allowed to use straightfor- 
wardly our Monte-Carlo simulation. The 68.3% confidence 
intervals plotted in Fig. || indicate that the PP-TSVD 
method can lead to an inferred width around 0.03i?© 
(which is the value obtained from LOWL data) for ini- 
tial widths up to O.O8i?0. Therefore the interpretation of 
the result obtained by this method is not easy. This may 
indicate that the method is better suited to the search of 
transition zones known a priori to be very thin (searching 
for a width lower than O.O5i?0 for example). Nevertheless, 
all the above discussions indicate 0.020 < w < 0.070i?o 
as a reasonable interval for the true width, deduced from 
PP-TSVD method. 

All these approaches are globally consistent but lead to 
a relatively large dispersion of the results. Therefore our 
final estimation of the width of the solar tachocline in the 
equatorial plane is: w = 0.05 ± 0.03i?©. This estimation 
is in very good agre emen t with the result obtained by 
Charbonneau et al. ( 1997 ) and rema ins compatible with 
the value given by Kosovichev ( 1996 ) (c.f. Table 2). 



Table 2. Comparison of our resu lts with previous forward 
analysis. Charbonneau et al. (1997) and our work are for the 
same LOWL dataset (2/26/94-2/25/96) whereas Kosovichev 
(|1996|) has used the 1986-90 BBSO datasets. 



6. Conclusions 

This work presents an analysis of the determination of the 
characteristics of the tachocline at the equator by three 
different inverse methods. They are applied to the inver- 
sion of the splittings of the sectoral modes estimated as 
the sum of the three first odd coefficients of the expan- 









w/Rq 


This work 


0.695 ± 


0.005 


0.05 ±0.03 


Charbonneau et al. 


0.704 ± 


0.003 


0.050 ±0.012 


Kosovichev 


0.692 ± 


0.005 


0.09 ±0.04 



sion of the splittings in orthogonal polynomials defined 
by Schou et al. ( [1994 ). Two different choices of regulariza- 
tion parameters, the GCV and L-curve criteria, have been 
compared. Tests with artificial rotation laws have shown 
that in all cases the GCV criterion is less sensitive to the 
error level than the L-curve one and gives better results 
with low bias and dispersions in the range 0.03 — 0.011-R© 
of searched widths. This choice of the GCV criterion is 



in agreement with Barett ( 1993 ) and Thompson (1992) in 
another context. Hansen ( 1992j ) has shown that the GCV 
criterion is less adapted to highly correlated errors than 
the L-curve one. Our work may indicate in turn that we 
can neglect, as it has been done, the unknown correlation 
in LOWL data. 

Concerning the thickness of the tachocline, it appears 
that the MTSVD and PP-TSVD inversions give respec- 
tively an upper and lower estimate while the Tikhonov 
method corrected by deconvolution gives the most reli- 
able determination. We have estimated the systematic ef- 
fect of the latitudinal dependence of the rotation in the 
convection zone on the determination of the thickness of 
the tachocline and the rotation in the radiative interior. 
We have shown how the performance of the methods will 
be improved by lowering the level of noise in the data. 

The methods have been applied to the LOWL two 
years dataset leading to an estimation of the position f c = 
0.695 ± 0.005i? Q and the thickness w = 0.05 ± O.O3i? of 
the equatorial tachocline. In addition, we have obtained an 
estimation of the equatorial rotation £Iq below the convec- 
tion zone and above 0.4i? Q such that: 427.5 < flo < 434.5 
nHz and Cli from the top of the convection zone up to 
0.8-R© such that Qi = 459. 0± 1.5 nHz. Assuming that the 
rotation in the radiative interior is independent of latitude, 
this leads to a ratio fio/^i between 0.93 and 0.95 which 
is intermediate between the two theoretical predictions. 

Our results for the location and thickness of the equa- 
torial tachocline are in agr eemen t with the forward anal- 
ysis of Charbonneau et al. ( 1997| ) and with those of Basu 
applied on BBSO and GONG datasets (Basu 1997 ) using a 
different parameterization of the tachocline. The forward 
analysis can be viewed as non-linear least-squares meth- 
ods (least-squares methods because of the use of the % 2 
criterion and non linear because of the models used for the 
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rotation profile) but using only a very few number of pa- 
rameters (Charbonneau et al. (1997) use six parameters, 
Basu ( L997[ ) three and Kosovichcv (1996) only two). This 
kind of methods depend thus strongly on our knowledge 
of the global rotation profile which can be reached only 
by inversion techniques. In particular, in the above-cited 
works the latitudinal dependence of the rotation is fixed 
(as in 1.5D inversions). In this work, we have tried to in- 
vestigate the amount of informations about the tachocline 
that we can extract directly from the global inversions 
without a-priori knowledge (except for the regularization) 
on the rotation profile. There are less assumptions in this 
approach, and thus the tachocline parameters may be less 
constrained. The fact that the two approaches lead to sim- 
ilar results indicates in turn that the hypothesis used in 
the forward analysis are probably not too strong and are 
well adapted to the problem of inferring the tachocline 
from actual data. 

One of the interest of this work was our first attempt 
to use an inverse method with non-linear regularization in 
helioseismic case. The PP-TSVD method leads to a very 
large dispersion of the results for widths above 0.05-R© and 
then is difficult to interpret with actual data. Some efforts, 
in future work, should be useful to improve this kind of 
methods and the interpretation of their results taking into 
account their non-linearity. 
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Schou for providing the LOWL data and the anonymous ref- 
eree for constructive comments. 



where r p , p — l..N p are fixed break points distributed ac- 
cording to the density of turning points of modes (Corbard 
et al., 1997). The matrix R is then defined by: 



R = (Ri; 



I i = l,N 



Ri 



K n l{r)(fp{r)dr 



(A5) 



For all the inverse methods discussed in this work, the 
aim is to find a solution that is able to produce a good fit of 
the data in chi-square sense. Unfortunately, the solution of 
this problem is not unique and allows oscillatory solutions 
that arc not physically acceptable. So, we have to define 
a quantity that measures the smoothness of the solution 
and to insure that the final solution is sufficiently smooth 
to be acceptable. 

For any solution f2, we define the % 2 value by: 



x 2 m = \\p(Rn-w)\\ 2 2 



(A6) 



where P = diag(l/o~ n i) and we define two measures of the 
smoothness of the solution 17 by: 



Pi{f2) = \\LS2\ 



i = 1,2 



(A7) 



where the vector i- norms ||.||j are defined by ||ce||i = 
(Ep kpl*) 1 ^ an d L is a discrete approximation of the first 
derivative operator such that: 



Pi (X 



9fi(r) 



Or 



dr 



/3|cx 



/ an(r) 

\ dr 



dr 



(A8) 



A.l. Tikhonov solution 



Appendix A: Details of the three inverse methods 
used 



We discretize Eq. (||) by: 
W = Rfl 

where we have defined: 



(Al) 



W = (W-)«=i,w Wi = Au nU + 5v nlh i = (n, I), (A2) 

N being the number of modes (n, I) (N = 1102 for LOWL 
data) and 5v n u a normally distributed noise with a stan- 
dard deviation defined in Eq. (|J). We search the solution 
f2(r) as a piecewise linear function of the radius by setting: 



P =i 



fl — (uip) p= i 



N„ 



(A3) 



where <p p (r), p = 1, N p are piecewise straight lines (N p 
50 in this work) such that: 



Vp = l.JVp, 3 r p e [0.,1.] / n(r p ) 



(A4) 



The so called Tikhonov solution fl\ solves the problem: 
min ( X 2 (f2) + Xf3 2 2 (f2)), (A9) 

where A > is the continuous regularization parameter. 
In order to compare this method to the two other ones, it 
may be interesting to reformulate the problem as follow: 
For any A we can show that there exist a value cv(A) for 
which f2\ is the solution of the problem: 

min (3 2 (f2); S x = {fl / \\P{Rf2 - W)\\ 2 < a(X)} (A10) 

The computation of these solutions for different regular- 
ization parameters have been carried out by using a gen- 
eralized singular value decomposition of the pair (J?, L) 
as explained and discussed extensively in Christensen- 
Dalsgaard et al. ( 1993 ). 



A.2. MTSVD and PP-TSVD solutions 

These methods are based on the SVD of the N x N p (N > 
N p ) matrix R which can be written: 



R = ^2 U ^i V I 

i=l 



(All) 
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where r < N p is the rank of R. The singular vectors are 
orthonormal, ujxij = 



vjv 



j = 5ij for i,j — l,r, and the 



singular values Aj are such that: Ai > A2 > ... > A r > 0, 
A r+ i, .., Ajv. = 0. We then define the the TSVD of R as 



the matrix Rk built from Eq. (All) but neglecting the 
N p — k smallest singular values. 



Rk 



u,.\,r 



(A12) 



The integer k < r is called the truncation parameter. It 
acts as a regularization parameter by eliminating the os- 
cillatory behavior of the singular vectors associated with 



the smallest singular values. According to Eq. (A12), the 
rank of the matrix Rk is k < r and then the problem 
of minimizing the quantity \\P(RkS2 — W)||a has not an 
unique solution and we have to use our smoothness cri- 
teria to select a physically acceptable solution among the 
set of solutions defined by: 



S k = {f2 j \\P{R k n - W)\\ 2 = minimum} 



(A13) 



With these notations, the so-called MTSVD solution fl\ 
is defined by: 



^ = argmin/3 2 (fi) 



(A14) 



whereas the so-called PP-TSVD solution fl p k is denned by: 



fl\ = arg min Pi(fl) 



nes k 



(A15) 



The algorithms for computing these solutions are pre- 
sented in Hansen et al. ( |l992 ) and Hansen & Mosegaard 
( |l99(i| ) respectively. 

We just recall some important properties of the PP- 
TSVD solution: For any k < r the vector Lfl p k has at the 
most k—1 non zero elements. As the matrix L is a discrete 
approximation of the first derivative, this means that the 
solution vector fl p k consists on kj, < k constant blocks. 
From Eq. ( |A3| ) it follows that the inferred rotation f2(r) 
itself is obtained as a piecewise constant functions with a 
maximum of k pieces. The kb — 1 break points of this so- 
lution are selected by the procedure among the N p initials 
break points r p . Therefore this inversion is able to pro- 
duce a discontinuous solution without fixing a-priori the 
location of the discontinuity. Finally, we note that the so- 
lution J?£ obtained by this non-linear method can always 
be computed by applying a matrix H, to the data but 
this matrix is also a function of the data i.e. H = H(W). 
Thus we have: 



n p = h(w)w. 



(Ai6) 
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